notebook09 - Statistical Analysis

rm(list=ls(all=TRUE))
load("~/Post-doc/Data/Total Merged Data File (July 24 2023).RData")
str(CC.TotalData)

Data for mean velocity, angular velocity, turn and heading angles each will be analyzed by a mixed model ANCOVA with chemical concentration, and light as fixed (treatment) effects, flow velocity as the co-variate and replicate as a random (blocking affect) with individual paths as the unit of analysis.

##To run ANCOVA in R load the following packages:
library(car)
library(compute.es)
library(effects)
library(ggplot2)
library(multcomp)
library(pastecs)
library(WRS2)

##If you are using lm or aov make sure that you set the contrasts using the "contrasts" function before doing either aov or lm. R uses non-orthogonal contrasts by default which can mess everything up in an ANCOVA. If you want to set orthogonal contrasts use:

contrasts(dataname$factorvariable)=contr.poly()# of levels, i.e. 3) ##example
contrasts(CC.TotalData$Chlorophyll)=contr.poly(7)
##then run your model as

model.1=aov(dv~covariate+factorvariable, data=dataname) ## example
mod1=aov(vel.flow~Flow.rate+Chlorophyll, data=CC.TotalData)
##To view the model use:

Anova(model.1, type="III") ## example
Anova(mod1, type="III")
##Make sure you use capital "A" Anova here and not anova. This will give results using type III SS.

summary.lm(model.1) ##will give another summary and includes the R-sq. output. Example
summary.lm(mod1)

posth=glht(model.1, linfct=mcp(factorvariable="Tukey"))  ##gives the post-hoc Tukey analysis. Example
posth=glht(mod1, linfct=mcp(Chlorophyll="Tukey"))
summary(posth) ##shows the output in a nice format.

##If you want to test for homogeneity of regression slopes you can also include an interaction term for the IV and covariate. That would be:

model=aov(dv~covariate+IV+covariate:IV, data=dataname) ## example

mod2=aov(vel.flow~Flow.rate+Chlorophyll+Flow.rate:Chlorophyll, data=CC.TotalData) ## not right?
##If the interaction term is significant then you do not have homogeneity.

Stats on each level

##Looking at other stats we could run
library(nlme)
library(lme4)
library(effects)


## example
d <- data.frame(
    Y = rnorm(48),  ### e.g. velocity
    subject = factor(rep(1:12, 4)),  ### e.g. D_V_T
    A = factor(rep(1:2, each=24)),   #### e.g. Light
    B = factor(rep(rep(1:2, each=12), 2)))   ### e.g. Flow.rate

summary(aov(Y ~ A*B + Error(subject/(A*B)), data=d))  # Standard repeated measures ANOVA

mode1 <- lme(Y~A*B, random = ~1|subject, data=d)  ## nlme

lmer(Y ~ A*B + (1|subject) + (1|A:subject) + (1|B:subject), data=d)  ### lme4



#########    our data using nlme package     #########

mode1 <- lme(pitch.perfect~Chlorophyll*Light, random= ~1|D_V_T, data=CC.TotalData)

mode2 <- lme(pitch.perfect~Light, random= ~1|D_V_T, data=CC.TotalData)

anova(mode1, mode2) #if it shows only minor improvement, no need to include missing factor

anova(mode2)


#########    our data using l4me package     #########

model1 <- lmer(vel.flow ~ Flow.rate*Chlorophyll + (1|D_V_T) + (1|Flow.rate:D_V_T) + (1|Chlorophyll:D_V_T), data=CC.TotalData)  ###

model2 <- lmer(vel.flow ~ Flow.rate + (1|D_V_T) + (1|Flow.rate:D_V_T), data=CC.TotalData)

anova(model1, model2) #if it shows only minor improvement, no need to include missing factor

anova(model1)

######### plotting effect of the models    ##########
plot(effect("turn.angle*Flow.rate",se=TRUE, confidence.level=.95, mode2))


###### interaction plot with Lattice extra #########
library(latticeExtra)
xyplot(CC.TotalData$vel.flow~CC.TotalData$turn.angle|CC.TotalData$Light*CC.TotalData$Flow.rate,  main="Scatterplots by Light and Flow Rate", 
   ylab="Vel.flow", xlab="Turn Angle") +
   layer(panel.ablineq(lm(y ~ x), r.sq = TRUE,
   at = 0.75, adj=1:5), style = 4)


############## interaction plots ###########################

interaction.plot(CC.TotalData$vel.flow, ## x axis
                 CC.TotalData$Flow.rate, ## factor for lines
                 CC.TotalData$turn.angle,  ## y axis
                 fun = mean,
                 type = c("b"), legend = TRUE,
                 trace.label = "Flow Rate",
                 ylim = c(-1,180),
                 main = "Lights On",
                 xlab = "vel.flow",
                 ylab = "turn angle",
                 col = CC.TotalData$Flow.rate, pch = c(1:9),
                 axes = TRUE)
#library(tidyr)
#library(dplyr)
#CC.TotalData %>% mutate_all(~replace(., is.nan(.), 0))
CC.TotalData[is.nan(CC.TotalData)] <- 0
range(CC.TotalData$log.vel.flow)
##########################################################

rm(mode1)

rm(mode2)


mod1 <- lm(vel.flow~Flow.rate*Chlorophyll, data = CC.TotalData)
summary(mode2)

             Kernel Density Plots

Turn Angles first…

################################ Chlorophyll ########################
levels(CC.TotalData$Chlorophyll)

[1] "No Chlorophyll"     "Low Chlorophyll"    "Medium Chlorophyll"
[4] "High Chlorophyll"  

Now looking at velocity in relation to flow instead….

### Kernel density plots ################### Swimming velocity
##range(CC.TotalData$vel.flow)
## plot(log.vel.flow)
##log.vel.flow <- log10(CC.TotalData$vel.flow[CC.TotalData$vel.flow > 0])
#str(CC.TotalData)
#  vels <- (CC.TotalData$v[CC.TotalData$D_V_T==ind[i]])  ## vels <- log10(vels[vels>0])  ## from notebook01, could be useful to try this instead
################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)

[1] "0"   "0.6" "3"   "5.9" "8.9"

Now for headings…..

rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8)

Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd6' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd7' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd8' not found

Now for pitch….

For each level of Flow

f0 <- CC.TotalData[CC.TotalData$Flow.rate=="0"]

Error in `[.data.frame`(CC.TotalData, CC.TotalData$Flow.rate == "0") : 
  undefined columns selected

For each level of Chlorophyll

### Kernel density plots ################### Turn Angles
par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)

kd1 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 0")

kd2 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="0.6"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="3"])
lines(kd3, col='yellow', lwd=2)

kd4 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="5.9"])
lines(kd4, col='orange', lwd=2)

kd5 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="8.9"])
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)




##kd1 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.3")

kd2 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="0.6"])
plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.3")
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.6")

kd2 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="0.6"])
##plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.6")
lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 5.5")

##kd2 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



##kd1 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 6")

kd2 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="0.6"])
plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 6")
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 7")

##kd2 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



##kd1 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 13")

kd2 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="0.6"])
plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 13")
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 19")

##kd2 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)


###########################################################################

#############################  Light  ######################################

###########################################################################

par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

kd1 <- density(ChlA0$turn.angle[ChlA0$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 0")

kd2 <- density(ChlA0$turn.angle[ChlA0$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA4.3$turn.angle[ChlA4.3$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 4.3")

##kd2 <- density(ChlA4.3$turn.angle[ChlA4.3$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA4.6$turn.angle[ChlA4.6$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 4.6")

##kd2 <- density(ChlA4.6$turn.angle[ChlA4.6$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA5.5$turn.angle[ChlA5.5$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 5.5")

##kd2 <- density(ChlA5.5$turn.angle[ChlA5.5$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(ChlA6$turn.angle[ChlA6$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 6")

##kd2 <- density(ChlA6$turn.angle[ChlA6$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA7$turn.angle[ChlA7$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 7")

kd2 <- density(ChlA7$turn.angle[ChlA7$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(ChlA13$turn.angle[ChlA13$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 13")

##kd2 <- density(ChlA13$turn.angle[ChlA13$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA19$turn.angle[ChlA19$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 19")

kd2 <- density(ChlA19$turn.angle[ChlA19$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

For each level of Flow (now for velocity)

### Kernel density plots ################### Velocity
par(mar= c(4,4,2,1))
par(mfrow=c(1,2))

################################ Chlorophyll ########################

kd1 <- density(f0$vel.flow[f0$Chlorophyll=="0"])
plot(kd1, xlim = c(-2,2), ylim = c(0, 15), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Chlorophyll, Flow = 0")

kd2 <- density(f0$vel.flow[f0$Chlorophyll=="4.3"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(f0$vel.flow[f0$Chlorophyll=="4.6"])
lines(kd3, col='yellow', lwd=2)

kd8 <- density(f0$vel.flow[f0$Chlorophyll=="5.5"])
lines(kd8, col='orange', lwd=2)

kd4 <- density(f0$vel.flow[f0$Chlorophyll=="6.1"])
lines(kd4, col='red', lwd=2)

kd5 <- density(f0$vel.flow[f0$Chlorophyll=="7.6"])
lines(kd5, col='dark red', lwd=2)

kd6 <- density(f0$vel.flow[f0$Chlorophyll=="13.5"])
lines(kd6, col='purple', lwd=2)

kd7 <- density(f0$vel.flow[f0$Chlorophyll=="19"])
lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chloropphyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)




kd1 <- density(f.6$vel.flow[f.6$Chlorophyll=="0"])
plot(kd1, xlim = c(-2,2), ylim = c(0, 15), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Chlorophyll, Flow = 0.6")

kd2 <- density(f.6$vel.flow[f.6$Chlorophyll=="4.3"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(f.6$vel.flow[f.6$Chlorophyll=="4.6"])
lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f.6$vel.flow[f.6$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

kd4 <- density(f.6$vel.flow[f.6$Chlorophyll=="6.1"])
lines(kd4, col='red', lwd=2)

##kd5 <- density(f.6$vel.flow[f.6$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

kd6 <- density(f.6$vel.flow[f.6$Chlorophyll=="13.5"])
lines(kd6, col='purple', lwd=2)

##kd7 <- density(f.6$vel.flow[f.6$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)



kd1 <- density(f3$vel.flow[f3$Chlorophyll=="0"])
plot(kd1, xlim = c(-15,15), ylim = c(0, 2.5), col='blue', lwd=2, main = "Velocity (bot logged, in relation to flow) Vs Chlorophyll, Flow = 3")

##kd2 <- density(f3$vel.flow[f3$Chlorophyll=="4.3"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(f3$vel.flow[f3$Chlorophyll=="4.6"])
##lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f3$vel.flow[f3$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

##kd4 <- density(f3$vel.flow[f3$Chlorophyll=="6.1"])
##lines(kd4, col='red', lwd=2)

##kd5 <- density(f3$vel.flow[f3$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

##kd6 <- density(f3$vel.flow[f3$Chlorophyll=="13.5"])
##lines(kd6, col='purple', lwd=2)

##kd7 <- density(f3$vel.flow[f3$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)



kd1 <- density(f5$vel.flow[f5$Chlorophyll=="0"])
plot(kd1, xlim = c(-15,15), ylim = c(0, 2.5), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Chlorophyll, Flow = 5.9")

##kd2 <- density(f5$vel.flow[f5$Chlorophyll=="4.3"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(f5$vel.flow[f5$Chlorophyll=="4.6"])
##lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f5$vel.flow[f5$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

##kd4 <- density(f5$vel.flow[f5$Chlorophyll=="6.1"])
##lines(kd4, col='red', lwd=2)

##kd5 <- density(f5$vel.flow[f5$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

##kd6 <- density(f5$vel.flow[f5$Chlorophyll=="13.5"])
##lines(kd6, col='purple', lwd=2)

##kd7 <- density(f5$vel.flow[f5$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)



kd1 <- density(f8$vel.flow[f8$Chlorophyll=="0"])
plot(kd1, xlim = c(-15,15), ylim = c(0, 0.25), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Chlorophyll, Flow = 8.9")

##kd2 <- density(f8$vel.flow[f8$Chlorophyll=="4.3"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(f8$vel.flow[f8$Chlorophyll=="4.6"])
##lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f8$vel.flow[f8$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

##kd4 <- density(f8$vel.flow[f8$Chlorophyll=="6.1"])
##lines(kd4, col='red', lwd=2)

##kd5 <- density(f8$vel.flow[f8$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

##kd6 <- density(f8$vel.flow[f8$Chlorophyll=="13.5"])
##lines(kd6, col='purple', lwd=2)

##kd7 <- density(f8$vel.flow[f8$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)



##########################################################################

################################# Guano ##################################

###########################################################################
par(mar= c(4,4,2,1))
par(mfrow=c(3,2))
levels(CC.TotalData$Guano)

kd1 <- density(f0$vel.flow[f0$Guano=="Present"])
plot(kd1, xlim = c(-10,10), ylim = c(0, 25), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Guano, Flow = 0")

kd2 <- density(f0$vel.flow[f0$Guano=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)

##kd1 <- density(f.6$vel.flow[f.6$Guano=="Present"])
##plot(kd1, xlim = c(-7,7), ylim = c(0, 0.12), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Guano, Flow = 0.6")

kd2 <- density(f.6$vel.flow[f.6$Guano=="Absent"])
plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Guano, Flow = 0.6")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)

##kd1 <- density(f3$vel.flow[f3$Guano=="Present"])
##plot(kd1, xlim = c(-7,7), ylim = c(0, 0.12), col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Guano, Flow = 3")

kd2 <- density(f3$vel.flow[f3$Guano=="Absent"])
plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Guano, Flow = 3")
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)

##kd1 <- density(f5$vel.flow[f5$Guano=="Present"])
##plot(kd1, xlim = c(-7,7), ylim = c(0, 0.12), col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Guano, Flow = 5.9")

kd2 <- density(f5$vel.flow[f5$Guano=="Absent"])
plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Guano, Flow = 5.9")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)

##kd1 <- density(f8$vel.flow[f8$Guano=="Present"])
##plot(kd1, xlim = c(-100,100), ylim = c(0, 0.04), col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Guano, Flow = 8.9")

kd2 <- density(f8$vel.flow[f8$Guano=="Absent"])
plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Guano, Flow = 8.9")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)


##########################################################################

################################# Light ##################################

###########################################################################
par(mar= c(4,4,2,1))
par(mfrow=c(3,2))
levels(CC.TotalData$Light)

kd1 <- density(f0$vel.flow[f0$Light=="Present"])
plot(kd1, xlim = c(-10,10), ylim = c(0, 25), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 0")

kd2 <- density(f0$vel.flow[f0$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(f.6$vel.flow[f.6$Light=="Present"])
plot(kd1, xlim = c(-10,10), ylim = c(0, 025), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 0.6")

kd2 <- density(f.6$vel.flow[f.6$Light=="Absent"])
##plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 0.6")
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(f3$vel.flow[f3$Light=="Present"])
plot(kd1, xlim = c(-10,10), ylim = c(0, 25), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 3")

##kd2 <- density(f3$vel.flow[f3$Light=="Absent"])
##plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 3")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(f5$vel.flow[f5$Light=="Present"])
plot(kd1, xlim = c(-10,10), ylim = c(0, 25), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Guano, Flow = 5.9")

##kd2 <- density(f5$vel.flow[f5$Light=="Absent"])
##plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 5.9")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(f8$vel.flow[f8$Light=="Present"])
plot(kd1, xlim = c(-10,10), ylim = c(0, 25), col='blue', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 8.9")

kd2 <- density(f8$vel.flow[f8$Light=="Absent"])
##plot(kd2, xlim = c(-10,10), ylim = c(0, 25), col='red', lwd=2, main = "Velocity (not logged, in relation to flow) Vs Light, Flow = 8.9")
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

For each level of Chlorophyll (velocity)

### Kernel density plots ################### velocity
par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

################################ Flow Rate ########################
kd1 <- density(log10(ChlA0$vel.flow[ChlA0$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 0")

kd2 <- density(log10(ChlA0$vel.flow[ChlA0$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA0$vel.flow[ChlA0$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(ChlA0$vel.flow[ChlA0$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA0$vel.flow[ChlA0$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)


kd1 <- density(log10(ChlA4.3$vel.flow[ChlA4.3$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 4.3")

kd2 <- density(log10(ChlA4.3$vel.flow[ChlA4.3$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA4.3$vel.flow[ChlA4.3$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(ChlA4.3$vel.flow[ChlA4.3$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA4.3$vel.flow[ChlA4.3$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(log10(ChlA4.6$vel.flow[ChlA4.6$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 4.6")

kd2 <- density(log10(ChlA4.6$vel.flow[ChlA4.6$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA4.6$vel.flow[ChlA4.6$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(ChlA4.6$vel.flow[ChlA4.6$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA4.6$vel.flow[ChlA4.6$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(log10(ChlA5.5$vel.flow[ChlA5.5$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 5.5")

kd2 <- density(log10(ChlA5.5$vel.flow[ChlA5.5$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA45.5$vel.flow[ChlA5.5$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(Chl5.5$vel.flow[ChlA5.5$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA5.5$vel.flow[ChlA5.5$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(log10(ChlA6$vel.flow[ChlA6$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 6")

kd2 <- density(log10(ChlA6$vel.flow[ChlA6$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA6$vel.flow[ChlA6$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(ChlA6$vel.flow[ChlA6$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA6$vel.flow[ChlA6$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(log10(ChlA7$vel.flow[ChlA7$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 7")

kd2 <- density(log10(ChlA7$vel.flow[ChlA7$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA7$vel.flow[ChlA7$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(ChlA7$vel.flow[ChlA7$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA7$vel.flow[ChlA7$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(log10(ChlA13$vel.flow[ChlA13$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 13")

kd2 <- density(log10(ChlA13$vel.flow[ChlA13$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA13$vel.flow[ChlA13$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(ChlA13$vel.flow[ChlA13$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA13$vel.flow[ChlA13$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(log10(ChlA19$vel.flow[ChlA19$Flow.rate=="0"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Flow Rate, Chlorophyll = 19")

kd2 <- density(log10(ChlA19$vel.flow[ChlA19$Flow.rate=="0.6"]))
lines(kd2, col='light blue', lwd=2)

kd3 <- density(log10(ChlA19$vel.flow[ChlA19$Flow.rate=="3"]))
lines(kd3, col='yellow', lwd=2)

kd4 <- density(log10(ChlA19$vel.flow[ChlA19$Flow.rate=="5.9"]))
lines(kd4, col='orange', lwd=2)

kd5 <- density(log10(ChlA19$vel.flow[ChlA19$Flow.rate=="8.9"]))
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)




#############################  Light  ######################################
par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

kd1 <- density(log10(ChlA0$vel.flow[ChlA0$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 0")

kd2 <- density(log10(ChlA0$vel.flow[ChlA0$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(log10(ChlA4.3$vel.flow[ChlA4.3$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 4.3")

kd2 <- density(log10(ChlA4.3$vel.flow[ChlA4.3$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(log10(ChlA4.6$vel.flow[ChlA4.6$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 4.6")

kd2 <- density(log10(ChlA4.6$vel.flow[ChlA4.6$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(log10(ChlA5.5$vel.flow[ChlA5.5$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 5.5")

kd2 <- density(log10(ChlA5.5$vel.flow[ChlA5.5$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(log10(ChlA6$vel.flow[ChlA6$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 6")

kd2 <- density(log10(ChlA6$vel.flow[ChlA6$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(log10(ChlA7$vel.flow[ChlA7$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 7")

kd2 <- density(log10(ChlA7$vel.flow[ChlA7$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(log10(ChlA13$vel.flow[ChlA13$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 13")

kd2 <- density(log10(ChlA13$vel.flow[ChlA13$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(log10(ChlA19$vel.flow[ChlA19$Light=="Present"]))
plot(kd1, col='blue', lwd=2, main = "Velocity (logged, in relation to flow) Vs Light, Chlorophyll = 19")

kd2 <- density(log10(ChlA19$vel.flow[ChlA19$Light=="Absent"]))
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


range(CC.TotalData$turn.angle)

Now for each level of Flow (for headings)

### Kernel density plots ################### Headings
par(mar= c(4,4,2,1))
par(mfrow=c(1,2))

f0$heading <- atan2(f0$dx, f0$dy)
range(f0$heading)
f0$heading.pi <- f0$heading/(2*pi)*360
range(CC.TotalData$heading)  ## gives -3.14 to 3.14
range(f0$heading.pi)  ## gives -180 to 180  


f.6$heading <- atan2(f.6$dx, f.6$dy)
range(f.6$heading)
f.6$heading.pi <- f.6$heading/(2*pi)*360
range(CC.TotalData$heading)  ## gives -3.14 to 3.14
range(f.6$heading.pi)  ## gives -180 to 180  


f3$heading <- atan2(f3$dx, f3$dy)
range(f3$heading)
f3$heading.pi <- f3$heading/(2*pi)*360
range(CC.TotalData$heading)  ## gives -3.14 to 3.14
range(f3$heading.pi)  ## gives -180 to 180  


f5$heading <- atan2(f5$dx, f5$dy)
range(f5$heading)
f5$heading.pi <- f5$heading/(2*pi)*360
range(CC.TotalData$heading)  ## gives -3.14 to 3.14
range(f5$heading.pi)  ## gives -180 to 180  


f8$heading <- atan2(f8$dx, f8$dy)
range(f8$heading)
f8$heading.pi <- f8$heading/(2*pi)*360
range(CC.TotalData$heading)  ## gives -3.14 to 3.14
range(f8$heading.pi)  ## gives -180 to 180  


################################ Chlorophyll ########################
levels(CC.TotalData$Chlorophyll)

par(mar= c(4,4,2,1))
par(mfrow=c(1,2))

kd1 <- density(f0$heading.pi[f0$Chlorophyll=="0"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Flow = 0", xlab = "Heading (degrees)")

##kd2 <- density(f0$heading.pi[f0$Chlorophyll=="4.3"])
##lines(kd2, col='light blue', lwd=2)

kd3 <- density(f0$heading.pi[f0$Chlorophyll=="4.6"])
lines(kd3, col='yellow', lwd=2)

kd8 <- density(f0$heading.pi[f0$Chlorophyll=="5.5"])
lines(kd8, col='orange', lwd=2)

##kd4 <- density(f0$heading.pi[f0$Chlorophyll=="6.1"])
##lines(kd4, col='red', lwd=2)

kd5 <- density(f0$heading.pi[f0$Chlorophyll=="7.6"])
lines(kd5, col='dark red', lwd=2)

##kd6 <- density(f0$heading.pi[f0$Chlorophyll=="13.5"])
##lines(kd6, col='purple', lwd=2)

kd7 <- density(f0$heading.pi[f0$Chlorophyll=="19"])
lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0",  "Chlorophyll 4.6", "Chlorophyll 5.5",  
                  "Chlorophyll 7.6",  'Chlorophyll 19'),
       col = c("blue",  "yellow", "orange",  "dark red", "black"),
       lty = 1)



kd1 <- density(f.6$heading.pi[f.6$Chlorophyll=="0"])
plot(kd1,xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Flow = 0.6", xlab = "Heading (degrees)")

kd2 <- density(f.6$heading.pi[f.6$Chlorophyll=="4.3"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(f.6$heading.pi[f.6$Chlorophyll=="4.6"])
lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f.6$heading.pi[f.6$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

kd4 <- density(f.6$heading.pi[f.6$Chlorophyll=="6.1"])
lines(kd4, col='red', lwd=2)

##kd5 <- density(f.6$heading.pi[f.6$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

kd6 <- density(f.6$heading.pi[f.6$Chlorophyll=="13.5"])
lines(kd6, col='purple', lwd=2)

##kd7 <- density(f.6$heading.pi[f.6$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 6.1", 
              "Chlorophyll 13.5"),
       col = c("blue", "light blue", "yellow",  "red", "purple"),
       lty = 1)


kd1 <- density(f3$heading.pi[f3$Chlorophyll=="0"])
plot(kd1,xlim = c(-180,180), ylim = c(0, 0.03), col='blue', lwd=2, main = "Headings Vs Chlorophyll, Flow = 3")

##kd2 <- density(f3$heading.pi[f3$Chlorophyll=="4.3"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(f3$heading.pi[f3$Chlorophyll=="4.6"])
##lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f3$heading.pi[f3$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

##kd4 <- density(f3$heading.pi[f3$Chlorophyll=="6.1"])
##lines(kd4, col='red', lwd=2)

##kd5 <- density(f3$heading.pi[f3$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

##kd6 <- density(f3$heading.pi[f3$Chlorophyll=="13.5"])
##lines(kd6, col='purple', lwd=2)

##kd7 <- density(f3$heading.pi[f3$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlrophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)



kd1 <- density(f5$heading.pi[f5$Chlorophyll=="0"])
plot(kd1,xlim = c(-180,180), ylim = c(0, 0.03), col='blue', lwd=2, main = "Headings Vs Chlorophyll, Flow = 5.9")

##kd2 <- density(f5$heading.pi[f5$Chlorophyll=="4.3"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(f5$heading.pi[f5$Chlorophyll=="4.6"])
##lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f5$heading.pi[f5$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

##kd4 <- density(f5$heading.pi[f5$Chlorophyll=="6.1"])
##lines(kd4, col='red', lwd=2)

##kd5 <- density(f5$heading.pi[f5$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

##kd6 <- density(f5$heading.pi[f5$Chlorophyll=="13.5"])
##lines(kd6, col='purple', lwd=2)

##kd7 <- density(f5$heading.pi[f5$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)

kd1 <- density(f8$heading.pi[f8$Chlorophyll=="0"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.03), col='blue', lwd=2, main = "Headings Vs Chlorophyll, Flow = 8.9")

##kd2 <- density(f8$heading.pi[f8$Chlorophyll=="4.3"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(f8$heading.pi[f8$Chlorophyll=="4.6"])
##lines(kd3, col='yellow', lwd=2)

##kd8 <- density(f8$heading.pi[f8$Chlorophyll=="5.5"])
##lines(kd8, col='orange', lwd=2)

##kd4 <- density(f8$heading.pi[f8$Chlorophyll=="6.1"])
##lines(kd4, col='red', lwd=2)

##kd5 <- density(f8$heading.pi[f8$Chlorophyll=="7.6"])
##lines(kd5, col='dark red', lwd=2)

##kd6 <- density(f8$heading.pi[f8$Chlorophyll=="13.5"])
##lines(kd6, col='purple', lwd=2)

##kd7 <- density(f8$heading.pi[f8$Chlorophyll=="19"])
##lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)


#####################################################################

################################# Guano ##################################

######################################################################

par(mar= c(4,4,2,1))
par(mfrow=c(3,2))
levels(CC.TotalData$Guano)

kd1 <- density(f0$heading.pi[f0$Guano=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.03), col='blue', lwd=2, main = "Headings Vs Guano, Flow = 0")

kd2 <- density(f0$heading.pi[f0$Guano=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)



##kd1 <- density(f.6$heading.pi[f.6$Guano=="Present"])
##plot(kd1, xlim = c(-180,180), ylim = c(0, 0.01), col='blue', lwd=2, main = "Headings Vs Guano, Flow = 0.6")

kd2 <- density(f.6$heading.pi[f.6$Guano=="Absent"])
plot(kd2, xlim = c(-180,180), ylim = c(0, 0.03), col='red', lwd=2, main = "Headings Vs Guano, Flow = 0.6")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)



##kd1 <- density(f3$heading.pi[f3$Guano=="Present"])
##plot(kd1, xlim = c(-180,180), ylim = c(0, 0.01), col='blue', lwd=2, main = "Headings Vs Guano, Flow = 3")

kd2 <- density(f3$heading.pi[f3$Guano=="Absent"])
plot(kd2, xlim = c(-180,180), ylim = c(0, 0.03), col='red', lwd=2, main = "Headings Vs Guano, Flow = 3")
##lines(kd2, col='red', lwd=2)

legend("topright",                              # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)



##kd1 <- density(f5$heading.pi[f5$Guano=="Present"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.03), col='blue', lwd=2, main = "Headings Vs Guano, Flow = 5.9")

kd2 <- density(f5$heading.pi[f5$Guano=="Absent"])
plot(kd2, xlim = c(-180,180), ylim = c(0, 0.03), col='red', lwd=2, main = "Headings Vs Guano, Flow = 5.9")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)



##kd1 <- density(f8$heading.pi[f8$Guano=="Present"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.03), col='blue', lwd=2, main = "Headings Vs Guano, Flow = 8.9")

kd2 <- density(f8$heading.pi[f8$Guano=="Absent"])
plot(kd2, xlim = c(-180,180), ylim = c(0, 0.03), col='red', lwd=2, main = "Headings Vs Guano, Flow = 8.9")
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)


################################################################

################################ Light ########################

#################################################################


par(mar= c(4,4,2,1))
par(mfrow=c(3,2))
levels(CC.TotalData$Light)

kd1 <- density(f0$heading.pi[f0$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.015), col='blue', lwd=2, main = "Headings Vs Light, Flow = 0")

kd2 <- density(f0$heading.pi[f0$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(f.6$heading.pi[f.6$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.015), col='blue', lwd=2, main = "Headings Vs Light, Flow = 0.6")

kd2 <- density(f.6$heading.pi[f.6$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(f3$heading.pi[f3$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.015), col='blue', lwd=2, main = "Headings Vs Light, Flow = 3")

##kd2 <- density(f3$heading.pi[f3$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(f5$heading.pi[f5$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.015), col='blue', lwd=2, main = "Headings Vs Light, Flow = 5.9")

##kd2 <- density(f5$heading.pi[f5$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(f8$heading.pi[f8$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.015), col='blue', lwd=2, main = "Headings Vs Light, Flow = 8.9")

kd2 <- density(f8$heading.pi[f8$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

Now for each level of Chlorophyll (for headings)

### Kernel density plots ################### Headings
par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)

par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

kd1 <- density(ChlA0$heading.pi[ChlA0$Flow.rate=="0"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.012), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 0")

kd2 <- density(ChlA0$heading.pi[ChlA0$Flow.rate=="0.6"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(ChlA0$heading.pi[ChlA0$Flow.rate=="3"])
lines(kd3, col='yellow', lwd=2)

kd4 <- density(ChlA0$heading.pi[ChlA0$Flow.rate=="5.9"])
lines(kd4, col='orange', lwd=2)

kd5 <- density(ChlA0$heading.pi[ChlA0$Flow.rate=="8.9"])
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)




##kd1 <- density(ChlA4.3$heading.pi[ChlA4.3$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 4.3")

kd2 <- density(ChlA4.3$heading.pi[ChlA4.3$Flow.rate=="0.6"])
plot(kd2, xlim = c(-180,180), ylim = c(0, 0.012), col='light blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 4.3")
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA4.3$heading.pi[ChlA4.3$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA4.3$heading.pi[ChlA4.3$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA4.3$heading.pi[ChlA4.3$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA4.6$heading.pi[ChlA4.6$Flow.rate=="0"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.012), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 4.6")

kd2 <- density(ChlA4.6$heading.pi[ChlA4.6$Flow.rate=="0.6"])
##plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 4.6")
lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA4.6$heading.pi[ChlA4.6$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA4.6$heading.pi[ChlA4.6$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA4.6$heading.pi[ChlA4.6$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA5.5$heading.pi[ChlA5.5$Flow.rate=="0"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.012), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 5.5")

##kd2 <- density(ChlA5.5$heading.pi[ChlA5.5$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA5.5$heading.pi[ChlA5.5$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA5.5$heading.pi[ChlA5.5$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA5.5$heading.pi[ChlA5.5$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



##kd1 <- density(ChlA6$heading.pi[ChlA6$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 6")

kd2 <- density(ChlA6$heading.pi[ChlA6$Flow.rate=="0.6"])
plot(kd2, xlim = c(-180,180), ylim = c(0, 0.012), col='light blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 6")
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA6$heading.pi[ChlA6$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA6$heading.pi[ChlA6$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA6$heading.pi[ChlA6$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA7$heading.pi[ChlA7$Flow.rate=="0"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.012), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 7")

##kd2 <- density(ChlA7$heading.pi[ChlA7$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA7$heading.pi[ChlA7$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA7$heading.pi[ChlA7$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA7$heading.pi[ChlA7$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



##kd1 <- density(ChlA13$heading.pi[ChlA13$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 13")

kd2 <- density(ChlA13$heading.pi[ChlA13$Flow.rate=="0.6"])
plot(kd2, xlim = c(-180,180), ylim = c(0, 0.012), col='light blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 13")
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA13$heading.pi[ChlA13$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA13$heading.pi[ChlA13$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA13$heading.pi[ChlA13$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA19$heading.pi[ChlA19$Flow.rate=="0"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.012), col='blue', lwd=2, main = "Heading Vs Flow Rate, Chlorophyll = 19")

##kd2 <- density(ChlA19$heading.pi[ChlA19$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA19$heading.pi[ChlA19$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA19$heading.pi[ChlA19$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA19$heading.pi[ChlA19$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)


###########################################################################

#############################  Light  ######################################

###########################################################################

par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

kd1 <- density(ChlA0$heading.pi[ChlA0$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 0")

kd2 <- density(ChlA0$heading.pi[ChlA0$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA4.3$heading.pi[ChlA4.3$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 4.3")

##kd2 <- density(ChlA4.3$heading.pi[ChlA4.3$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA4.6$heading.pi[ChlA4.6$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 4.6")

##kd2 <- density(ChlA4.6$heading.pi[ChlA4.6$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA5.5$heading.pi[ChlA5.5$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 5.5")

##kd2 <- density(ChlA5.5$heading.pi[ChlA5.5$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(ChlA6$heading.pi[ChlA6$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 6")

##kd2 <- density(ChlA6$heading.pi[ChlA6$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA7$heading.pi[ChlA7$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 7")

kd2 <- density(ChlA7$heading.pi[ChlA7$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(ChlA13$heading.pi[ChlA13$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 13")

##kd2 <- density(ChlA13$heading.pi[ChlA13$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA19$heading.pi[ChlA19$Light=="Present"])
plot(kd1, xlim = c(-180,180), ylim = c(0, 0.02), col='blue', lwd=2, main = "Headings Vs Light, Chlorophyll = 19")

kd2 <- density(ChlA19$heading.pi[ChlA19$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)